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Abstract 

Quantitative photoacoustic tomography (QPAT) is a recent hybrid imaging modal¬ 
ity that couples optical tomography with ultrasound imaging to achieve high resolution 
imaging of optical properties of scattering media. Image reconstruction in QPAT is 
usually a two-step process. In the first step, the initial pressure field inside the medium, 
generated by the photoacoustic effect, is reconstructed using measured acoustic data. 

In the second step, this initial ultrasound pressure field datum is used to reconstruct 
optical properties of the medium. We propose in this work a one-step inversion al¬ 
gorithm for image reconstruction in QPAT that reconstructs the optical absorption 
coefficient directly from measured acoustic data. The algorithm can be used to recover 
simultaneously the absorption coefficient and the ultrasound speed of the medium from 
multiple acoustic data sets, with appropriate a priori bounds on the unknowns. We 
demonstrate, through numerical simulations based on synthetic data, the feasibility of 
the proposed reconstruction method. 

Key words. Photoacoustic tomography, hybrid inverse problems, image reconstruction, one-step 
reconstruction, numerical optimization. 


1 Introduction 

Photoacoustic tomography (PAT) is a recent multi-physics biomedical imaging modality 
that aims at achieving simultaneously high resolution and high contrast in imaging by cou¬ 
pling the high-resolution ultrasound imaging modality with the high-contrast diffuse optical 
tomography (DOT) modality. In PAT, near infra-red (NIR) photons are sent into an op¬ 
tically absorbing and scattering medium, for instance a piece of biological tissue, Q C 
(d > 2), where they diffuse. The density of the photons, denoted by w(x), solves the following 
diffusion equation [5, 9] 

—V ■ D(x)Vii(x) + cr(x)«(x) = 0, in Q , , 

u = g(x), on <9Q 
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where D and cr are respectively the diffusion and absorption coefficients of the medium, and 
g is the incoming photon source. The medium absorbs a portion of the incoming photons and 
heats up due to the absorbed energy. The heating then results in thermal expansion of the 
tissue and the expansion generates a pressure field. This process is called the photoacoustic 
effect. The pressure field generated by the photoacoustic effect can be written as [9, 21] 

H(x) = T(x)cr(x)u(x) (2) 


where er(x)w(x) is the total energy absorbed locally at x G 12 and T is the Griineisen 
coefficient that describes the efficiency of the photoacoustic effect. This initial pressure field 
H then propagates in the form of ultrasound with sound speed c(x) [9, 21, 57] 


1 d 2 p 
c 2 (x) dt 2 


— — A p 

) t 2 y 

= 0, 

in M + x il 

p{ 0,x) 

= H(X), 

in 12 

'x 

CD 

&H 1 +■ 

\CZ 

= 0, 

in f2 

p(t,x) 

= 0, 

on M + x <9f2 


( 3 ) 


assuming that the medium has no acoustic attenuation effect. We refer interested readers 
to [7, 9, 2 ] for the derivation and justification of the models (1) and (3). 

The time-dependent acoustic signals that arrive on the surface of the medium, |(o ,T)xdci, 
are then measured with acoustic devices for a sufficiently long time T. From the knowledge 
of these acoustic measurements, one is interested in reconstructing the diffusion and absorp¬ 
tion properties of the medium; see [6, 12, 16, 34, 36, 39, 46, 53, 63, 64] for overviews of 
photoacoustic tomography. 

Image reconstructions in PAT are usually done in a two-step process. In the first step, 
one uses the boundary acoustic signal to reconstruct the initial pressure field H in (2), 
inside the medium. This step has been extensively studied in the past decade under various 
circumstances; see for instance [1, 2, 4, 13, 14, 20, 24, 25, 27, 33, 35, 37, 44, 46, 48, 57, 60] and 
references therein. In the second step, usually called the quantitative step, one reconstructs 
the diffusion coefficient D (x), the absorption coefficient cr(x), and whenever possible, the 
Griineisen coefficient T(x) using the internal datum i/(x). This step has recently received 
great attention as well; see for instance [3, 7, 9, 10, 15, 17, 22, 38, 41, 43, 47, 50, 52, 55, 68] 
and references therein. Uniqueness and stability results for the reconstruction procedures 
have been established in different settings. 

The above two-step process works perfectly fine in the setting where acoustic measure¬ 
ment can be performed everywhere on the boundary <9f2 and the ultrasound speed c(x) is 
known. In this case, the initial pressure field F7(x) can be reconstructed relatively stably in 
the first step of PAT reconstructions. In the case where acoustic data can only be measured 
on part of the boundary, i.e. the so called limited-view setting, the reconstruction of H in 
the first step can be fairly unstable (especially in the part of the domain from where acoustic 
signals do not travel easily to the measurement locations) [14, 2 ]. In the case where the 
sound speed c(x) is also unknown, which is often the case in real-world applications, the 
first step reconstruction is also problematic. In this case, one has to reconstruct simulta¬ 
neously the wave speed c and the initial pressure field H. The reconstruction has been 
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shown recently by Stefanov and Uhlmann [58], based on a slightly different formulation of 
the acoustic problem, to be very unstable. 


In practical applications, we almost always measure acoustic data generated from multi¬ 
ple optical illuminations. The two-step process, however, does not take advantage of these 
multiple data sets. The reason is that if we change the illumination source g, the initial 
pressure field H will also be changed. Therefore, every time we add a new measurement, we 
introduce a new unknown H in the first step of the PAT reconstruction. 

The ultimate objectives of the PAT reconstructions are the coefficients (c, D, cr, F). These 
coefficients do not change when the illumination is changed. Therefore adding more data 
should not introduce more unknowns in the reconstruction. In other words, instead of 
reconstructing the unknowns {c,H) and then the coefficients (D,a, T), we should recon¬ 
struct directly from the measured data the coefficients (c, D , cr, T), without the intermediate 
variable H. This way, one can potentially use data sets from multiple illuminations to help 
stabilize the reconstruction, even in the case of partial measurements (for each illumination). 


The aim of the current work is exactly to develop one-step reconstruction strategies that 
recover the optical coefficients and the wave speed directly from the measured acoustic data. 
For simplicity of presentation, we focus only on the absorption coefficient cr and the wave 
speed assuming that the diffusion coefficient D and the Griineisen coefficient T are known. 
In this case, we have to invert the nonlinear map A(c, cyg) defined through the following 
relation: 


dp 

7^|(o,T)xan = A(c, cr; g). 


(4) 


To the best of our knowledge, there is no uniqueness result on this inverse problem of 
reconstructing simultaneously the ultrasound speed and the absorption coefficient, besides 
the special case in [33]. The result in [58] indicates that the reconstruction would be unstable 
when only one data set, i.e. data collected from only one optical illumination, is used. Our 
main goal here is to show numerically that using multiple data sets allows us to obtain 
fairly stable reconstructions, with appropriate a priori bounds on the unknown coefficients. 
Rigorous mathematical study of the stability of our reconstruction approach will be a future 
work. 


To set up the problem, we assume in the rest of the paper that (A-i) the domain is 
bounded with smooth boundary dfl, (A-ii) the boundary condition g is the restriction of a 
C 4 function on dFl, (A-iii) the coefficients c(x) G C 4 (f2), cr(x) G C 4 (fl), -D(x) G C' 3 (k2), and 
T(x) G C 4 (k2), and (A-iv) (c, a, D, T) G A a x Ag x A$ x Ag with 

A a = {/(x) : 0 < a < /(x) < a < oo, Vx G fl}, (5) 

a and a being two constants. With the assumptions (A-i)-(A-iv), we conclude from standard 
elliptic theory [19, 23] that the diffusion equation (1) admits a unique solution u G C 4 ( fl). 
This implies that Fan G C 4 {Fl). Therefore the wave equation (3) admits a unique solution 
p G C 4 ((0,T) x f2) following the theory in [18, ‘ ]. The datum (4) is then well-defined and 
can be viewed as a map A : C 4 (f2) x C 4 (fl) i-G "H 5//2 ((0,T) x dFl). 

The rest of the paper is structured as follows. We first present a one-step inversion algo¬ 
rithm in the linearized setting in Section 2. We then implement in Section 3 the algorithm 
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for reconstructions in the fully nonlinear setting. We show some numerical reconstructions 
based on synthetic data in Section 4 to demonstrate the feasibility of the algorithm. Con¬ 
cluding remarks are offered in Section 5. 


2 Inversion with Born approximation 


We denote by J the total number of illumination sources available, and ^-|(o,T)x 9 n the 
measured datum on the boundary for illumination g 3 (1 < j < J), in the time interval (0, T). 
Our objective is thus to numerically invert the following nonlinear system to reconstruct 
(c,a): 

^ l~in l(o,T)x9f2 \ ( A(c, a; gi) \ 


dp 

dn 


V *rl(0 ,T)xdn J 


= A(c,<r;g) 


( 6 ) 


A (c,a-,gj) 

In this section, we develop a one-step algorithm for the inverse problem in the linearized 
setting where we intend to reconstruct perturbations to the coefficients around known back¬ 
grounds. This is useful in some real-world applications where variations of the ultrasound 
speed and the absorption coefficient are relatively small (for instance, it is well-known that 
ultrasound speed in biological tissues is very similar to that in water with about 15% vari¬ 
ations from tissue to tissue [67]). 


2.1 The Born approximation 


We denote by c 0 (x) and cr 0 (x) respectively the known background ultrasound speed and 
absorption coefficient. We assume that the true coefficients are of the forms 

c(x) = c 0 (x) + c(x), and cr(x) = cr 0 (x) + d(x), (7) 

where the perturbations c and a satisfy ||c/co||L°°(n) •C 1 and ||dt/cr 0 1|•C 1 respectively. 
The perturbations in the coefficients lead to perturbations in the solutions to the diffusion 
equation and the wave equation as follows: 


u j = u °j + Uj , and pj (t , x) = p ° (t, x) + pj (t, x) 
where the background photon densities are determined as solutions to 

—V ■ DX7u°j(-x) + a 0 (x)u°(x) = 0, in 

9ji 


11 — on dil 


Uj = 


( 8 ) 


( 9 ) 


and the background acoustic pressure fields p^ solve the acoustic wave equations: 


d 2 p° j 


Cq(x) dt 2 


-A p°j = 0, 


in (0, T) x VL 


P?(0,x) = r<T 0 ri°(x), in Q 
.0 


dp\ 

-qP (°> x ) = 

P°(t,x) = 0, 


( 10 ) 


in 12 

on (0, T ) x <912 






We check that the perturbations Uj and pj solve respectively, after neglecting higher-order 
terms, the diffusion equation: 


—V ■ UVuj{x.) + <7 0 (x)iij(x) = —cfu° p in f] 

Uj = 0, on <9f] 


and the wave equation: 

1 d 2 pj 


Apj = 


2 < d 2 p°j 


-c x 


Cq(x) dt 2 3 Cq dt 2 

Pj( 0,x) = r(du° + (T 0 Uj), inf] 


in (0, T) x f] 


A(0,x) = 0, 


in f] 

on (0, T ) x <9f] 


Pj(t,x) = 0, 

Let us now introduce the linear operators A^(co,<7o) and A£(co,oo) through: 

-|(0,T)x<9Q = A;)(c 0 , CTo)d, 

un 

with p c - and pj respectively the solutions to 


dpj 


- 7 ^\(o,T)xdn = A J c (c 0 ,u 0 )c, and ^ 


1 d 2 p‘ 


and 


£T - A P C J 

2 d 2 pj 

= 3 0+2 C ( X )> 

Cq dt 2 

in (0, T) x f] 

P% 0,x) 

= o, 

in f] 

'x 

o' 

= 0, 

in f] 


= 0, 

on (0, T ) x <9f] 

b 

< 

1 

h 1 

= 0, 

in (0, T) x f] 


1 d 2 p c j 
Cq(x) dt 2 

pj(0,x) = r(du° + a 0 Uj), inf] 

dp° 

J (0,x) = 0, 


in f] 

on (0, T ) x <9f] 


(ii) 


( 12 ) 


(13) 


(14) 


(15) 


dt 

Pj(t,x ) = 0, 

We can then write the perturbation of the datum as 

dp ■ 

-^(t,x) |( 0 ,T)xen = A^(c 0 ,u 0 )c +A^(co,u 0 )d. 

We then collect data for all J sources to get a system of equations for the unknowns: 

c \ <9p 


A CjCT (co, (To) 


a 


dn 


with 


A-c,cr(Ccb ° o ) 


1 A)(c 0 ,(t 0 ) 

A ^(coWo) \ 
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(16) 


(17) 


(18) 
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This is the Born approximation of the original nonlinear problem (6). The approximation 
can be justified as stated in the following lemma. 

Lemma 2.1. Let 12, D, T and gj satisfy the assumptions in (A-i)-(A-iv). Then the datum 
generated from gj, viewed as the map: 

dp ■ 

A(c,<7; 9j ): (c ’ <t) in l(0 ' T|x8n (19) 

C 4 (f2) x C\Q) xdn) 

is Frechet differentiable at any (c 0 ,<r 0 ) G C 4 (f2) x C 4 (f2) that satisfies the assumption in 
(A-iv). The derivative at (c 0 ,cr 0 ) in the direction (c, <r) G C 4 (f2) x C 4 (f2) (such that c 0 + c 
and (Jo + a satisfy (A-iv)) is (A^(c 0 , cr 0 )c, A£(co, <7o)<j) as defined in (13). 

The lemma can be proven using standard arguments such as those in [11, 18, 29, 49]. 
We provide a sketch of the proof in the Appendix. Note that even though A(c, a; gf) is 
well-defined as a map A(c, cr; gf) : C 4 (f2) x C 4 (12) i —y "H 5//2 ((0, T) x <912), we can only prove its 
differentiability as a map A(c, cr\gf) : C 4 (12) x C 4 (12) i-> "H i / 2 ((0,T) x <912). 

We now need to invert (16) (or (17) when multiple data sets are available) to reconstruct 
the perturbation (c,<r). It is well-known that with a single measurement ^(f, x)[( 0 ,T)xan, 
one could reconstruct one of c and Hj = r(<T?j( J + a 0 uf) (thus a since Hj uniquely determines 
a [8, 9]) assuming that the other is known [29, 59, 65, 66]. In a slightly different setting, 
Stefanov and Uhlmann [ ] showed that even if a single measurement ^-(f, x)|( 0 ,T)xan is 

enough to reconstruct the pair (c, Hf (and thus (c, a)) uniquely, the reconstruction would 
be extremely unstable [58, Theorem 1], Our hope here is that, without introducing new 
unknowns (since we reconstruct <j directly, not Hf, by using multiple data sets we can 
improve the stability of the reconstruction, assuming that uniqueness can be achieved. 


2.2 Reconstruction based on Born approximation 


To invert the linear system (17), we use the technique of Landweber iteration [32], The 
iteration takes the form: 


{ Qc+1 
V &k +1 


(I - tA* ><t (c o, <to)A C)CT (c 0 , ct 0 )) 


Ck 


+ T -^*c,a( C 0: °~o) 


dp 

dn’ 


k > 0, 


( 20 ) 


with a reasonable given initial guess. The parameter r, 0 < r < 2/E 2 with E being the 
largest singular value of A C)f7 (c 0 , <t 0 ), is a positive algorithmic parameter that we select by 
trial and error since we do not have good estimates on the singular values of A CjO .(c 0 , cr 0 ). 
The components of the adjoint operator A* a (co,<Jo) 


Ac,cr( C 0> °b) 


Ac*( C 0 ; CTo) 
Aa*( C 0; CTo) 


A c*(Co,0-o) \ 
A/*(co,<Jo) ) ' 


( 21 ) 


are given as follows. The adjoint operator A/*(c 0 ,<Jo), in the sense of 

(Ac(c 0 ,<J 0 )c, Uj)L 2 (gd,T)xdn) = {c,Ai*{c 0 ,a 0 )y j } L 2 {n) , G T 2 ((0,T) x dfl), 
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is given as 


K*( c o,vo)yj = 


r ' T 2 d 2 p° j 


dt 2 


(t, x)dt 


( 22 ) 


with qj the solution to the following adjoint wave equation: 

1 d 2 q 3 


Cq(x) dt 2 


- Aqj = 0, 
dj(T,x) = 0, 
f (T,x) = 0, 


in (0, T) xQ 
in hi 
in hi 

on (0, T) x dQ 


(23) 


«(*, x ) = -Vi, 

The adjoint operator A£*(co, cr 0 ), in the sense of 

(A^(c 0 , a 0 )cr, Zj) L 2 (( 0 T ^ xd ty = (a,Ai*(c 0 ,a 0 )z j ) L 2 {n) , \/ Zj e T 2 ((0,T) x dfl), 
is given as 

K'*s = _ (^ r ^(° )X ) + Vj) u °ji 

where qj is now the solution to (23) with y 3 replaced by Zj and v 3 is the solution to the 
adjoint diffusion equation: 


(24) 


1 do 

V ■ D'Vvjix) + cr 0 (x)vj(x) = — 2 r ^o^-(0,x), in ff 

Cn Oi 


(25) 


= 0, 


on <9ff 


Note that even though we split the components for c and er in the datum in the form 
of (16) for the convenience of presentation, we do not need to solve two wave equations, (14) 
and (15), to compute the datum. Instead we need only to solve one wave equation, i.e (12). 
Therefore, in each iteration of the Landweber algorithm, we need to solve J forward wave 


equations and J forward diffusion equations to evaluate f*. = A C>t7 (c 0 , cr 0 ) 


Ck 


and then J 


adjoint wave equations and J adjoint diffusion equations to evaluate A* CT (co, cro)ffc. The last 
term in the iteration does not change during the iteration, so it needs only to be computed 
once before the iteration starts. 


3 One-step nonlinear reconstruction 


To solve the full nonlinear inverse problem, we take an optimal control approach. We look 
for the solution to the inverse problem as 


Cmin 

^"min 


= argmin F(c, a), 

(c,a)eA a xAp 


(26) 


with the objective functional F given as: 


<9p # 


F (c,a) = -||A(c,a;g) - ^ ||[L2( (0 ,Dxan)]-'> 
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(27) 




where is the collection of measured acoustic data. 

We implemented the Levenberg-Marquardt method [31, 45] to solve the minimization 
problem. The method is characterized with the following iteration: 


f Qc+i 

\ &k +1 


C-k 

®k 


(K,A c <c^t)^cA c ^ a t) + m) l K.A Ct *•■ > o. 


(28) 


where /i*. is an algorithm parameter, is the residual at step k\ 


Zfc = A(c fc , a k ] g) 


a P # 

<9n 


and A Cj(T (cfc, cr*,) is the Frechet derivative of A at (c*,, cr*,). In our implementation, we take the 
Levenberg-Marquardt parameter fi k as a small constant, although we are aware that there 
are principles in the literature to guide the selection of this parameter in a more “optimal” 
way; see for instance [31, 45]. 

The Levenberg-Marquardt algorithm (28) requires the inverse of the operator (A^Ac^-l- 
fi k I) at each iteration. In our implementation, we do not form the operator (which in discrete 
case is a matrix) explicitly and then invert it, since this would require large computer 
memory to store the matrix. Instead, we use a matrix-free approach to save memory in 
the following way. For any function (which in discrete case is a vector) y, to compute 
2 = (K Ac, a + Hk I) V, we solve the following equation: 


(A^A^ + fi k I)z = y. 


(29) 


This is a symmetric positive definite problem which we solve with a standard conjugate 
gradient method [ ]. The conjugate gradient method does not require the explicit form 

of the operator {A* ca A C)(y + /i/J), but only its action on given functions (vectors in discrete 
case). For any given z, we solve J forward wave equations and J forward diffusion equations 
to evaluate f = A CiCr (cfc, cqjz and then J adjoint wave equations and J adjoint diffusion 
equations to evaluate A* CT (cfc, cp^f. 

To impose the bound constraints on the coefficients, that is c G A a (i.e. a < c(x) < a) 
and a € Ap (i.e. f3 < cr(x) < /3), we implement the algorithm for the new variables i'(x) 
and ?/(x) that are related respectively to c(x) and cr(x) through the relations: 


. . a + a a — a . . 

c(x) =-1-tanhn(x), 


/3 + /3 P — /3 

cr(x) = =-1-=tanh? 7 (x), 


(30) 


The Frechet derivatives of A with respect to the new variables at (r>o, ?/o) hi the direction 
( v,fj ) can be computed using the chain rule straightforwardly as 


A„(t> 0 )S = A e (c(ti 0 ))[ Q 11 li ' I,ix)i 1 . A,(i) 0 )ij = A„(ct(i) 0 ))[—-^sech 2 i; 0 (x)i)]. 

(31) 

Note that we have chosen different bounds for c (i.e a and a) and a (i.e. (3 and (3) since in 
practice the two functions have different ranges of values. 








4 Numerical experiments 


We now present some numerical experiments to check the performance of the reconstruction 
algorithm. To simplify the presentation, we non-dimensionalize the problem so that all the 
numbers presented below have no dimensions. We consider the problem in the square domain 
0 = [0 2] x [0 2] with constant diffusion coefficient D = 0.02 and Griineisen coefficient T = 1. 

In our implementation, we discretize the forward and adjoint diffusion equations, for 
instance (1) and (25), with a first-order finite element method. We discretize the forward 
and adjoint wave equations, for instance (3) and (23), with a standard second-order finite 
difference scheme on a uniform grid. The finite element discretization of the diffusion equa¬ 
tions is performed on a triangle mesh that shares the same nodes as the uniform grid for 
the wave equation. This way we do not need to interpolate between the two types of grids 
when using quantities from the diffusion solution in the wave equations or vice versa. 

We will perform numerical reconstructions in both the nonlinear and the linearized set¬ 
tings. To generate synthetic data for the nonlinear inversions, we solve the diffusion equa¬ 
tion (1) and the wave equation (3) with the true absorption coefficient and ultrasound speed, 
and then compute the data using (4). When adding multiplicative random noise to the da¬ 
tum we perform the transformation ^ x (1 + ^ • rand) where rand is a uniform 

random variable with mean 0 and variance 1 (and thus range [— \/3, V^3])- We use n to 
measure the level of noise in the data. The values of k will be given in the simulations we 
present below. For the linearized inversions, we construct synthetic data directly using the 
linearized model (16) for a given perturbation of the coefficients, (c, a). This means that 
the data for the linearized inversions are exact. These data do not contain information on 
the accuracy of the linearized problem as an approximation to the true nonlinear problem. 
The linearized inversion simulations we show below will only provide information about the 
invertability and stability of the linearized inverse problem. 

We measure the quality of the reconstruction with the maximal relative error. For 
parameter p, the error is defined as ||(p r — p*)/p* ||x,°° where p 4 is the true coefficient while p r 
is the reconstructed coefficient. 

Numerical experiment 1. The first numerical experiment is devoted to the reconstruc¬ 
tion of the absorption coefficient assuming that the ultrasound speed is known. The true 
absorption coefficient is taken as 



0.15, x e [0.5 1.5] x [0.5 1.5] 
0.10, x G T2\[0.5 1.5] x [0.5 1.5]. 


(32) 


We performed nonlinear reconstructions with the Levenberg-Marquardt algorithm using 
three types of data: (i) noise-free data (n = 0.0), (ii) noisy data with k = 0.5, and (iii) noisy 
data with k = 1.0. The data are collected from eight different optical illuminations that 
are line sources supported on each side of the domain with spatially varying strengths such 
as those used in [8]. The results of the reconstructions are shown in Fig. 1. The maximal 
relative errors in the reconstructions are 0.15, 0.28 and 0.64 respectively. The quality of the 
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reconstructions is very similar to those published in the literature [3, 15, 17, 22, 38, 41, 43, 
47, 50, 52, 55, 68]. We observe in our numerical experiments that the reconstructions in this 
case are very robust to changes in initial guesses. Moreover, we impose very loose bounds 
on the absorption coefficient in the reconstructions: /3 = 0 < a < 1.0 = /3. Our experience 
is that this bound is not necessary at all when a is the only unknown to be reconstructed. 






Figure 1: Reconstructions of the absorption coefficient with known ultrasound speed. Left 
to right: true cr(x) given in (32), a reconstructed from noise-free data, a reconstructed from 
noisy data with n = 0.5, and a reconstructed from noisy data with k = 1.0. 


Numerical experiment 2. In the second numerical experiment, we reconstruct the ul¬ 
trasound speed c(x) assuming that the absorption coefficient is known. We consider the 
problem with the true ultrasound speed 

lx — 11 llP 

c(x) = 1.0 + 0.2 x exp(--— ’ ), X e Cl. (33) 

Z X u.o 

We again performed reconstructions with the three types of data (i)-(iii). The results are 
shown in Fig. 2. We imposed the bound a = 0.8 < c < 1.3 = a on the unknown. The max¬ 
imal relative errors in the reconstructions are respectively 0.16, 0.30 and 0.57. The initial 
guess for all the reconstructions shown is c = 0.9, although we observe in our numerical ex¬ 
periments that most constant initial guesses within the bounds lead to almost identical final 
reconstructions. The data used in the reconstructions are again collected from 8 different 
optical illuminations. Even though we have more blurring in the reconstructed images, the 
overall quality of the reconstructions is reasonable, considering that we do not need very 
strict bounds on the unknown to get these reconstructions. 


Numerical experiment 3. In this experiment, we reconstruct, under the Born approx¬ 
imation, the sound speed c and the absorption coefficient a as shown in the Erst column 
of Fig. 3. The background of the linearization is (co,<7o) = (1.0, 0.1). We again use data 
collected from eight different illuminations. Shown are the reconstructions from noise-free 
data (k = 0.0, column 11), the reconstructions from noisy data with /t = 0.5 (column III) 
and the reconstructions from noisy data with n = 1.0 (column IV). The maximal relative 
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1.051.1 1.15 1.051.11.151.2 1 1.051.11.15 



Figure 2: Reconstructions of the ultrasound speed with known absorption coefficient. Left to 
right: true c(x) given in (33), c reconstructed from noise-free data (k = 0.0), c reconstructed 
from noisy data with k = 0.5, and c reconstructed from noisy data with k = 1.0. 


errors in the reconstructions are (0.40,0.90), (0.80,1.20), and (1.05,2.02) respectively. We 
started all the Landweber iterations with initial condition (c, a) = (0,0). Let us empha- 





»«i 





0.08 0.1 0.12 



0.08 0.1 0.12 



0.08 0.1 0.12 



0.080.1 0.12 


Figure 3: Linearized reconstructions of the coefficients c (top) and a (bottom). Shown from 
left to right are: the true (c, a), the reconstructions from noise-free data, the reconstructions 
from noisy data with k = 0.5, the reconstructions from noisy data with k = 1.0. 

size again that the synthetic data used in this experiment are constructed directly from the 
linearized model (16). Therefore, the data are exact (besides the artificial noise added to 
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them) in the sense that the error in the approximation to the true nonlinear problem is 
completely neglected. What we are interested in studying is the stability of the linearized 
inverse problem, not the accuracy of the Born approximation. This is why we can consider 
fairly large perturbations to the background here. 


Numerical experiment 4. We repeat the reconstructions in the previous experiment (i.e. 
Numerical experiment 3) in the nonlinear setting with the Levenberg-Marquardt algorithm. 
The results are shown in Fig. 4. The maximal relative errors in the reconstructions are 
(0.42,0.69), (0.61,1.02), and (0.89,1.87) respectively for data with noise level k = 0.0, 
k = 0.5 and n = 1.0. In all the reconstructions, we use the initial guess of (c, a) = (0.9, 0.09) 
and the Levenberg-Marquardt algorithm is stopped after 80 iterations. We impose the 
bounds a = 0.7 < c < 1.3 = a, (3 = 0.07 < cr < 0.15 = f3. 
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Figure 4: Same as Fig. 3 except that the reconstructions are performed in the nonlinear 
setting using the Levenberg-Marquardt algorithm. 


Numerical experiment 5. In the last numerical experiment, we reconstruct the sound 
speed c and the absorption coefficient a as shown in the first column of Fig. 5. We again use 
data collected from eight different illuminations. Shown are the reconstructions from noise- 
free data (k = 0.0, column II), the reconstructions from noisy data with k = 0.5 (column III) 
and the reconstructions from noisy data with k = 1.0 (column IV). The maximal relative 
errors in the reconstructions are (0.97,0.56), (1.17,0.88), and (1.41,0.99) respectively. The 


12 









Levenberg-Marquardt algorithm is stopped at iteration 100 in all the reconstructions in 
this case and the initial guess for all the reconstructions is (c, cr) = (1.0,0.11). We impose 
stricter bounds on c in this case to have the algorithm converge to reasonable solutions: 
a = 0.85 < c < 1.2 = a, /3 = 0.1 < a < 0.2 = /3. Note that the bounds on the absorption 
coefficient coincide with the bounds on the true cr value. If the constant initial guess is 
beyond the bounds, the algorithm usually does not converge. 



Figure 5: Nonlinear reconstructions of the coefficients c (top) and cr (bottom). Shown from 
left to right are: the true (c, cr), the reconstructions from noise-free data (k = 0.0), the 
reconstructions from noisy data with n = 0.5, the reconstructions from noisy data with 

K = 1 . 0 . 


We observe from all the simulations presented in this section that the one-step recon¬ 
struction strategy performs reasonably well in either the linearized or the nonlinear case, 
with a priori bounds imposed on the unknowns. Tuning algorithmic parameters, such as 
the r in the Landweber iteration and the jik in the Levenberg-Marquardt algorithm, the 
maximal number of iterations allowed for the iterative algorithms, and the initial guesses 
etc, can certainly improve the quality of the reconstructions slightly as we observed in some 
of the simulations we performed. We did not pursue in that direction. The use of data from 
even more illuminations would certainly help at least to reduce the average noise level in the 
data. We did not pursue in that direction either. The simulations we present here provide 
us a rough idea about the quality of the reconstructions that we could get. 


13 






















5 Concluding remarks 


There have been several works in recent years on the simultaneous reconstruction of the 
ultrasound speed and the optical properties [30, 33, 40, 42, 51, 58, 61, 62, 69] in QPAT. Every 
method proposed attempts to follow the two-step philosophy, that is to first reconstruct 
ultrasound speed and the initial pressure field and then reconstruct the optical properties. In 
some cases, additional ultrasound measurements are taken to supplement the photoacoustic 
measurements [42]. We proposed here a reconstruction strategy that combines the two-step 
reconstruction process into a one-step process to reconstruct directly the sound speed and the 
optical properties without reconstructing the initial pressure field, which is an intermediate 
variable. 

The main advantage of our method is that it allows the use of data sets from multiple 
illuminations which can stabilize the reconstruction when the ultrasound speed is treated as 
a unknown. When the intermediate variable, the initial pressure field, is to be reconstructed 
as in [30, 33, 40, 54, 58, 61, 62, 69], it changes with illuminations. Therefore adding data 
from more illuminations simply adds more unknowns in the reconstruction process. While 
in our case, adding more data sets does not add more unknowns since (c, a) do not change 
with illuminations. 

The results in [8, 9, 29, 59, 65, 66] show that one can reconstruct uniquely and stably 
either the ultrasound speed or the optical absorption coefficient if the other one is known. 
We are not aware of any uniqueness result on the simultaneous reconstruction of both 
coefficients besides in special cases such as these in [ 53 , 40]. Numerical simulations with 
synthetic data in Section 4 show that the simultaneous reconstruction is in fact relatively 
stable when multiple data sets are used, with appropriate a priori bounds imposed on the 
unknowns. This does not contradict the theoretical results in [58] on the instability of the 
problem with a single measurement. We are currently conducting theoretical studies on the 
stability property of our method. Results will be reported elsewhere. 

Let us finish this paper with two additional remarks. First, there have been similar one- 
step algorithms in the literature [26, 56] that reconstruct directly optical properties from 
photoacoustic data. These algorithms are different from the one we proposed here since they 
all treat the ultrasound speed as known. Second, in our formulation, we have assumed that 
both the diffusion coefficient D and the Griineisen coefficient T are known. We can easily 
modify the algorithm to include D or T as a unknown in the reconstructions as well. Note 
that it has been proved [8] that one can not reconstruct simultaneously D, a and T with 
mo no-wavelength optical illuminations. Therefore, it is not desired to reconstruct (c, D, a, T) 
in the one-step algorithm. However, one can indeed attempt to reconstruct simultaneously 
c, a and T, for instance. 
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Appendix: Frechet differentiability of A(c, a; gj) 

We provide a brief proof of the differentiability of A(c, cr,gj) as stated in Lemma 2.1. To 
simplify the presentation, we will use the short notation u — and therefore u> 0 — 4, 

C C 0 

u — — We will show differentiability of A(cu, a; gj) with respect to (co, a). 

Proof of Lemma 2.1. We first prove the differentiability with respect to oj (and thus c). Let 
Pj(u> o + a), (To) and Pj(oJo, cr 0 ) = pj be the solution to the wave equation with coefficients 

(cuo + cD, (To) and (a;o,<To) respectively. De&ie pj = Pj(coo + c5, cr 0 ) ~ Pj(wo,co) an d Pj = 
Pj(uJo Tcu, cr 0 ) — Pj(oJo, do) — pf , where pf is the solution to (14) with replaced by 0. We 
then verify that pj solves 


and pj solves 


(o; 0 + u 


d 2 pj 


8t 2 


Apj 

Pj{ 0,x) 


pj(t,x) 


d 2 p°- 

in (0 ,T)xQ 

0, in O 

0, in O 

0, on (0, T ) x <90 


d 2 Pj A x 

~ Ar 


Pj( °>x) 



(0,x) 


Pj(t,x) 


in (o,T)xO 
0, in O 

0, in O 

0, on (0, T) x <90. 


(34) 


(35) 


With the assumptions in the lemma, we conclude from standard elliptic theory [19, - 
that the diffusion equation (9) admits a unique solution if- G C 4 (0). This implies that 
r a 0 Uj G C 4 (0). Therefore the wave equation (10) admits a unique solution pj = pj(u> 0 , <To) G 

C 4 ((0,T) x O) following theory in [18, 29], Moreover, G C 2 ((0,T) x O) and we have 
from (34) that p 3 G C 3 ((0,T) x O) and 


_ d 2 pfj d 2 pj 

\\Pj\\n 3 {(o,T)xn) < ^'i\\^-Q^-\\n 2 ({o, T )xa) IIh 2 ((o,t)xq)||^II« 4 (o) 

< C , i||p 3 || W 4(( 0)T ) xO )||d;|| W 4(Q) < C' 2 ||r(ToW°||^4(Q)||d;||- W 4( f2 p (36) 
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Similarly, we have from (35) that 


^ ~ d 2 p ■ ~ d 2 p • 

»j||w 2 ((0,T)xO) < C 'l||^-^/llw 1 ((0,T)xn) < Cl 11 


|w 1 ((o,T)xfi)|p||w 4 (n) 

< Ci||pj|| W 3(( 0iT ) xQ )||a;|| W 4(Q). (37) 


We now combine the trace theorem, (36) and (37) to get 

c)pj ~ o 2 

ll _ ^HIw 1 / 2 ((o,T)xan) < C 2 ||r^ W 4 (n) ||<h|| w 


(38) 


This shows the differentiability with respect to to. 

To prove its differentiability with respect to a, we observe that A(c, a; gj) = A(c, Hj(a ); g.,) 
is linear with respect to Hj. Therefore, A is differentiable with respect to Hj e C 4 (0) with 

~ ~ dpH 

the derivative at (cq,Hj) in the direction Hj given as Ah, (cq, H { -)H 3 = -^|(o,T)xan where 


pf solves 


1 d 2 pf 


Apf = 


dt 2 


dp 


0, in (0, T) x 

pf( 0,x) = Hj, in O 
ii 

(0,x) = 0, in Q 

0, on (0, T ) x dQ. 


(39) 


dt 


We now recall that Hj(a) : C 4 (Q) —> C 4 (Q) is Frechet differentiable with the derivative at 
a 0 in direction a given as H ]a (a a )a = Y^au® + a 0 Uj) where u 3 solves (11). The chain rule 
of differentiation then concludes that A is differentiable with respect to a at <7o and the 
derivative is A£(c 0 , cr 0 )a = A Hj (c 0 , cr 0 )H j(7 (a 0 )d. 

□ 
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